Density-functional theory for attraction between like-charged plates 
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We study the interactions between two negatively charged macroscopic surfaces confining positive 
counterions. A density-functional approach is introduced which, besides the usual mean-field inter- 
actions, takes into account the correlations in the positions of counterions. The excess free energy 
is derived in the framework of the Debye-Hiickel theory of the one-component plasma, with the ho- 
mogeneous density replaced by a weighted density. The minimization of the total free energy yields 
the density profile of the microions. The pressure is calculated and compared with the simulations 
and the results derived from integral equations theories. We find that the interaction between the 
two plates becomes attractive when their separation distance is sufficiently small and the surface 
charge density is larger than a threshold value. 



I. INTRODUCTION 

Solutions containing macromolecules are ubiquitous in 
the everyday life. From food colloids to the DNA, we 
are surrounded by these giant molecules which directly 
or indirectly govern every aspect of our lives. In many 
cases the macromolecules in solution posses a net charge. 
The electrostatic repulsion between the polyions is, of- 
ten, essential to stabilization of colloidal suspensions. In 
the biological realm the electrostatics is responsible for 
the condensation of the DNA and formation of actin bun- 
dles, while various physiological mechanisms depend on 
the electrostatic interactions between the proteins and 
the microions. In spite of their ubiquity, our understand- 
ing of polyelectrolyte solutions is far from complete. 

The effort to fathom the role of electrostatics as it 
applies to the colloidal suspensions goes back over half 
a century to the classic works of Derjaguin_and Lan- 
dauEl and of Verwey and Overbeek (DLVO).H These in 
turn wereubased on the pioneering studies of GouyH and 
Chapmaro of double layers in metal electrodes. Follow- 
ing these early contributions, a large effort has been de- 
voted to solve the Poisson-Boltzmann (PB) equation in 
various geometries. The mean-field treatment, based on 
the solution of the PB equation, suggests that the in- 
teraction between two equally charged macroions irhpi 
suspension containing counterions is always repulsive Jj'EI 
In recent years, however, this dogma began to be qHOS-i 
tioned based on simulations,!! analytical calculationsO^tll 
and experiments,E3"E3 which indicated that for small 
distances and large charge densities, two like-charged 
polyions might actually attract! 

The fundamental goal of this paper is to demonstrate 



that this attraction is linked to the correlations between 
the microions omitted in the mean-field theories, and to 
establish the conditions under which the attraction be- 
comes possible. We shall consider the interaction be- 
tween two infinite uniformly charged plates confining 
their own point-like counterions. The mean-field approx- 
imation for this system is obtained by solving the PB 
equation which, due to the planar symmetry, can be done 
analytically. Once the density profile is obtained, all the 
other thermodynamic quantities can be easily derived. 
Thus, it is not difficult to demonstrate that the pres- 
sure at the mean-field level, in units of energy, is simply 
the density of counterions at the mid-plane between the 
plates. Since this is always positive, no attraction is pos- 
sible within the mean-field theory. 

Realization that the correlations between the coun- 
terions can strongly modify the mean-field predictions 
goes back a number of years. One of the £rst ap- 
proaches proposed by Kjellander and MarceljaH was to 
include the correlations through the numerical solution of 
the Anisotropic Hypernetted Chain Equation (AHNC). 
These authors found that the force per unit area (pres- 
sure) can become negative in the presence of divalent 
counterions. Monte |£!arlo (MC) simulations performed 
by Guldbrand et alB also indicate that as the surface 
charge density is increased, the pressure decreases if 
the distance between the charged surfaces is sufficiently 
small. As in the case of the AHNC calculations, attrac- 
tion was found only in the presence of divalent counteri- 
ons. These authors, however, did not analyze the case of 
very high charge density and short distance between the 
plates. In addition, since in the above calculations it is 
difficult to separate the different physical contributions 
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to the pressure, the mechanism that drives the attraction 
remains unclear. 

A different theoretical approach which attempted to 
shed some light on the mechanisrp^f attraction was ad- 
vanced by Stevens and RobbinsO These authors pro- 
posed a density-functional theory similar to the one of- 
ten employed in studies of simple liquids. This approach 
introduces a grand-potential free energy, [p(r)], which 
is a functional of the non-uniform density of counteri- 
ons p(r). The equilibrium properties of the system are 
obtained through the minimization of the total free en- 
ergy. The practical problem with this method is that the 
exact form of the functional is not known. When the 
correlations between the microions are omitted, the min- 
imization of the grand potential, Qpb, becomes trivial 
and leads to the usual PB equation.B In order to account 
for the cpjjrelations between the counterions, Stevens and 
Robbi 




■pealed to the Local Density Approximation 
(LDA)Jl3^Eil Within this approach an additional contri- 
bution, /lda, is added to the mean- field expression, flpB- 
The expression for /lda adopted by Stevens and Robbins 
was obtained through the extrapolation of the MC data, 
for the homogeneous One-Component Plasma (0CP),E2I 
but with the homogeneous density replaced by an inho- 
mogeneous density profile. The minimization of the free- 
energy functional allowed them to determine the density 
profile, p{r), and the pressure, Pqcp- The LDA, how- 
ever, is not without its own problems. The major draw- 
back of this approach is that, for short distances and high 
charge densities, the LDA is unstable. The reason for the 
instability is due to the fact that as the density of counte- 
rions in the vicinity of the plates increases, the chemical 
potential decreases, what attracts more particles to the 
region. This, in turn, leads to an unphysical "chain reac- 
tion" where all the counterions condense onto the plates. 
Clearly, when the distance between the counterions be- 
comes smaller than some threshold vaUie_*™rr, the LDA 
ceases to be a reliable approximation.EJ'E^lEJ 

An improvement over the LDA is, thjC, sp , called, 
Weighted Density Approximation (WDA).E50B In this 
case, the excess free energy is taken to be a function of 
an average density, Pw{r) = J d^r' w{\r — r'\) p{r'), aver- 
aged over a region of radius s — Scom where theinterac- 
tions between the counterions are the strongest B3'cll The 
difficulty in the practical implementation of this scheme 
is the determination of a proper weight function. The 
simplest possible, form for w{\r — r'|), used by Stevens 
and Robbins jl^a was to assume that this function haa-a, 
long-range variation comparable to the wall separation.E3 
In this case, the weighted density Pwi^) is approximated 
by the homogeneous density independent of r. However, 
when the walls are not close, L > Scom the weighted func- 
tion is no longer uniform and the approximation adopted 
by Stevens and Robbins becomes unrealistic. 

A beautiful explanation of the attraction between like- 
charged plates lias been recently advanced by Rouzina 
and Bloomfield.O These authors present a picture of 
attraction as arising from the ground-state configura- 



tion of the counterions. Clearly at zero temperature 
the counterions will recondense onto the surface of the 
plates forming two intercalating Wigner crystals. The 
authors advance a hypothesis that even at finite temper- 
atures, relevant to the common experimental conditions, 
the attraction is still governed by the zero-temperature 
correlations. A somewhat different formulations based 
on field-theoretic methodology have also been proposed. 
In these approximations the attraction arises as a re- 
sult of carjCfilated fluctuations in the counterion charge 
densities-LitJa Although providing a nice qualitative ex- 
planation of the origin of the attraction, these simple 
theories fail to yield a quantitative agreement with the 
simulations. 

In this paper we propose a different form of the 
weighted-density approach, which rectifies the problems 
of the earlier theories while still remaining numerically 
tractable. The excess free energy and the weight func- 
tion, w{\r — r'\), are deriveiLfrom the Debye-Hiickel-Hole 
(DHH) theory of the 0CP.E3 The density profile is deter- 
mined by minimizing the free-energy density with respect 
to the local density. Once the density profile is obtained, 
the free energy of the system is calculated by inserting it 
into the expression for the free-energy functional. Given 
the free energy, all the thermodynamic properties of the 
system can be easily calculated. A careful analysis of 
the behavior of the pressure as a function of the charge 
density and the distance between the plates allows us to 
explore the nature and the origin of the attraction. 

The remainder of the paper is organized as follows. 
The model and the PB approximation for the density- 
functional approach are described in Sec. 0. The WD A 
is introduced and applied in Sec. pil|. Our results and 
conclusions are summarized in Sec. IIVI- 



II. THE POISSON-BOLTZMANN APPROACH 

We consider two large, charged, thin surfaces each of 
area A, separated by a distance L (see Fig. ^. The two 
plates with a negative surface charge density, —a, confine 
positive point-like monovalent counterions with charge e. 
The overall charge neutrality of the system is guaranteed 
by the constraint 
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dzp{z) ^ — , 
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where p{z) is the local number density of counterions and 
z is the Cartesian coordinate perpendicular to the plates. 
The space between the plates is assumed to be a dielectric 
continuum of constant e. 

In order to explore the thermodynamic properties of 
the system, we use a density-functional approach. The 
grand potential of the system is 



n [p] = T[p] - fiN 



(2) 
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where N is the total number of counterions, /i is their 
chemical potential and the functional is derived from 
the free-energy density of the homogeneous system, with 
the uniform density of counterions, pc — N/LA, replaced 



by the local density p{z). For dilute systems, the ionic 
correlations can be neglected and the grand-potential 
functional (per unit area) becomes 
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where the electrostatic potential, 

£ r — r'\ 



(4) 



due to the symmetry of the problem, depends only on 
the z coordinate. A is the de Broglie thermal wave- 
length of the counterions, (3 = l/ksT and q{z) = 
—a [S{z — L/2) + 6{z + L/2)] is the surface charge den- 
sity of the plates. The functional minimization of this 
expression. 
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ASp{z) 

produces the optimum density profile, 
p{z) = po cxp [-/3e0(z)] 



(5) 



(6) 



The electrostatic potential is obtained by solving the 
Poisson equation. 
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with the distribution of free ions given by Eq. 
find 
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where (/)o is the reference potential, which we will set to 
zero. Here A = 1/\/2ttXbPo and As — /3e^/e is the Bjer- 
rum length. Eq. (||) has to obey two boundary conditions, 
namely. 
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FIG. 1. Two infinite, negatively charged thin plates, with 
surface charge density —a separated by distance L. The coun- 
terions are confined to the region between the plates. The sol- 
vent is modeled as a uniform medium of dielectric constant 



The constant po is determined from the overall charge- 
neutrality condition, Eq. (pi). 
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S(z = 0) = , 

L\ Ana 



(10) 



From the first equation, the electric field vanishes at the 
mid-plane and, therefore, zq = 0. The second equation 
imposes the discontinuity of the electric field at both 
charged surfaces, leading to 
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The potential at a point z is, then, given by 
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(11) 
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with A the root of Eq. ( pi] ) . The optimum density profile 
derived from this potential. 
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cos2(z/A) 



(13) 



can now be substituted into the free-energy functional, 
allowing the calculation of the total free energy. The 
thermodynamic properties of the system can be deter- 
mined from a suitable differentiation of the total free en- 
ergy. For example, the force between the two plates is 
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given by the minus derivative of the free energy with re- 
spect to the separation L between the two surfaces. This 
differentiation leads to a particularly simple expression 
for the force per unit of area (or pressure), 



/3P = Po 



(14) 



We note that although it might be tempting to attribute 
this simple result to the contact theorem, this is not the 
case, since the conditions under which this theorem holds 
are violated in the present geometry; Eq. ( p^ is purely 
a mean- field result. 



III. THE WEIGHTED-DENSITY 
APPROXIMATION 

For dense systems, the correlations between the mi- 
croions become relevant. For instance, if a counterion is 
present at the position r, due to the electrostatic repul- 
sion, the probability that another counterion is located 
in its vicinity is drastically reduced. The correlations in 
the positions of the counterions reduce the mean-field es- 
timate of the electrostatic free energy. No exact method 
exists for calculating this excess contribution. The sim- 
plest approximation, the LDA, consists of adding to the 
Eq. (0) a local functional, 
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/lda = / dz p{z)fcorT[p{z)] 
' -L/2 



(15) 



where /corr [p{z)] is the correlational free energy per par- 
ticle. Within the LDA one normally uses the expression 
derived for the homogeneous system, in which the uni- 
form density pc = N/ LA is replaced by the local density 
profile p(t). Unfortunately, as was mentioned above, the 
LDA is unstable when the one-particle density p{r) is a 
rapidly varying function of the position. For example, 
for high surface charge densitieSy-ithe minimization of the 
grand potential has no solution.cS To circumvent this and 
related problems intrinsic to the LDA, TarazonaEj and 
Curtin and Ashcroftn3 proposed a WDA, in which the 
free-energy density, /lda, is replaced by 
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The fundamental difference between the LDA and the 
WDA is that the latter is assumed to depend not on the 
local density p{r), but on some average density within 
the neighborhood of the point r. 



;(r) = J d^r'i 



]r~r'\;p{r)] p{r') 
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This provides a control mechanism which prevents an un- 
physical, singular, buildup of concentration at one point. 
The grand potential is obtained by adding the excess free 
energy per area, given by Eq. (16), to the Eq. (^, 
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Minimization of this expression leads to the optimum 
particle number density. 



p{z) = po exp [-/3e(/)(z) - Ppc^{z)] 



(19) 



where the excess chemical potential derived from /wda, 
Eq. dH), is 
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and the normalization coefficient is 
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The electrostatic potential satisfies the Poisson equation, 
Eq. (H), with the charge density given by the Eq. (|is|). 
Integrating the Poisson equation over a rectangular shell 
of area A and width z, and appealing to the Gauss' theo- 
rem, an integro-differential equation for the electric field 
E{z) can be obtained. 
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(22) 

where £ = ejSXgE, a = crAg/e, z = z/Xb, L = L/Xg 
and /icx = f3pcx- The local density p{z), which en- 
ters in the calculation of the excess chemical potential, 
Eq. ( pO[ ) , can be obtained from the derivative of the elec- 
tric field, since V • E{r) = AT:ep{r)/e. The Eq. (H) 
explicitly fulfills the two boundary conditions: £{0) = Q 
and £{±L/2) = ±4Tra. 

The solution of this equation depends on the specific 
form of the excess free-energy density and the weight 
function w{\r — r'\). For the homogeneous OCP the 
electrostatic free energy can be easily obtained using the 
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DHH theory of Nordholmll This is a simple hnear the- 
ory based on the ideas of Debye and Hiickel. The elec- 
trostatic potential of the OCP is assumed to satisfy a 
linearized PB equation. As a correction for the lineariza- 
tion, Nordholm postulated the existence of an excluded- 
volume region of size Scom from which all other ions are 
excluded. The size of this region is such that the electro- 
static repulsion between two counterions is comparable 
to the thermal energy. Recent calculations using a gen- 



eralized Debye- Hiickel theory indicate that this exclusion 
region is responsible for the oscillations obseriisd in the 
structure factor of the OCP at high couplingsO Follow- 
ing Nordholm, we find 



Scorr = (1 + SAbKd)^^'"^ 



(23) 



where kd = y/^nXsPc is the inverse of the Debye length. 
The excess free energy per particle is calculated to be 
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where to = {1 + SXbkdY^^ ■ The correlational free en- 
ergy per particle for the WDA, /com which appears in 
(pO|), is obtained by replacing by Pwiz) in the expres- 
sion (11), that is, /coir [Pw{z)] = fopWrjfic ^ Pw{z)]. 

To obtain the weighted functionHEj we require that 
the second functional derivative of the free energy J- in 
the limit of homogeneous densities, 
(52/3jC- S^{r-r 
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(25) 



produces the direct correlation function C2{r) of the ho- 
mogeneous system. 
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Following Groot,E3 we find that a reasonable approxima- 
tion for the weight function is 



w{r) — w{\r\) 
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e(sc 



r) , (27) 



where Q{x) is the Heaviside step function. It is impor- 
tant to remember that the radius of the excluded- volume 
region, Scom is now a function of the position, since the 
average density pc, which appears in Eq. (p^), is replaced 
by p{z), the local density of counterions, see Eq. (|T^). 
Taking advantage of the planar symmetry of the system, 
the expression for the weighted density can be written 
explicitly as a one-dimensional quadrature. 
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where z< = max(— i/2, z — Scorr), Zy — min(L/2, z 
Scorr) and Scorr is a function of z through p{z). 



IV. RESULTS AND CONCLUSIONS 



be determined from the numerical iteration of Eq. ( p^ ) 
until convergence is obtained. The Helmholtz free energy, 
_F, associated with the optimum countcrion distribution 
( |l9| ) , is determined by substituting it into the free-energy 
functional T , 



Once /corr, Mox(2) and pw{z) are defined, the electric 
field and, consequently, the optimum density profile can 
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Using this expression, the pressure, for different dis- 
tances between the plates, L, and various charge densi- 
ties, cr, can be easily obtained through numerical differ- 
entiation, 



AdL 



(30) 



as shown in Fig. When the charge density is below 
a threshold value, a < ac, the dimensionless pressure, 
Xgj3P, is always positive and a monotonically decreasing 
function of L. Above the critical surface charge density 
the pressure exhibits a distinct minimum. In particular, 
we find that for sufficiently high surface charge densities 
the force between the two like-charged surfaces becomes 
negative, i.e. the two plates attract! 




FIG. 2. The reduced osmotic pressure as a function of 
the plate separation for various surface charge densities 
a = aX%/e: 1 (A), 5 (□) and 7 (O)- The solid line is the 
WDA and the dashed line is the PB approximation for the 
same values of a. 



In order to compare our results with other theoriesJS'EJ 
we assumed that the dielectric medium between the 
plates is water at room temperature and, consequently, 
that the Bjerrum length is As = 7. 14 A. The distance 
between the plates is fixed at 150 A and the inverse of 
the surface charge density, S = e/a, is varied from 40 A^ 
to 1000 A^. Our results, illustrated in Fig. |[ show that 
for small surface charge densities the pressure increases 
almost linearly with the inverse charge density, S. In 
this case, since Pcorr ^ PpB, the pressure is dominated 
by the PB behavior. However, when the charge density 



becomes large, the slope of Pwda increases due to the 
strong repulsion between the countcrions. 
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FIG. 3. Variation of /3P with E = e/cr for L = 150 A: from 
PB (dashed), WDA (solid) and AHNC (•) from Ref. [8]. 

We also compare gur calculations with the simulations 
of Guldbrand et alB In this case, the distance between 
the plates is fixed at 21 A and the surface charge densit; 
is varied from 0.01 C/m^ to 0.6 C/m^, as shown in Fig 
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FIG. 4. The osmotic pressure as a function of the surface 
charge density, when the distance between the plates is fixed 
at i = 21 A, from PB (dashed) and WDA (solid). The circles 
(•) are the data from Ref. [7]. 

When the density of counterions is small, Pwda does not 
differ significantly from Prb- As the surface charge den- 
sity is increased, the correlations among the counterions 
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become relevant and Pwda changes its slope and begins 
to decrease. Our results are in good agreement with the 
simulations, which also indicate that for a separation of 
21 A the pressure exhibits a region where, it decreases 
with increase in the surface charge density.Q 



ACKNOWLEDGMENTS 

We acknowledge the fruitful discussions with Marcelo 
Louzada-Cassou, Rudi Podgornik, Roland Kjellander 
and Stjepan Marcelja. One of us, Marcia Barbosa, is 
particularly grateful for the useful discussion with Mark 
O. Robbins. This work was supported in part by CNPq 
— Conselho National de Desenvolvimento Cientffico e 
Tecnologico and FINEP — Financiadora de Estudos e 
Projetos, Brazil. This research was also supported by the 
National Science Foundation under Grant No. PHY94- 
07194. 



^ B. V. Derjaguin and L. Landau, Acta Physicochimica 
(USSR) 14, 633 (1941). 

^ E. J. W. Verwey and J. Th. G. Overbeek, Theory of the Sta- 
bility of Lyophobic Colloids (Elsevier, Amsterdam) 1948. 

^ G. Gouy, J. Phys. 9, 457 (1910). 

* D. L. Chapman, Philos. Mag. 25, 475 (1913). 

^ J. N. Israelachvili, Intermolecular and Surface Forces, 2nd 
Edition (Academic Press, London) 1992. 

® S. A. Safran, Statistical Thermodynamics of Surfaces, In- 
terfaces and Membranes (Addison- Wesley, Reading Mass.) 
1994. 

L. Guldbrand, B. Jonsson, H. Wennerstrom and P. Linse, 
J. Chem. Phys. 80, 2221 (1984). 



* R. KjeUander and S. Marcelja, Chem. Phys. Lett. 112, 49 

(1984); J. Phys. Chem. 90, 1230 (1986). 
3 R. Podgornik, J. Chem. Phys. 91, 5840 (1989). 
M. Lozada-Cassou and E. Diaz-Herrera, J. Chem. Phys. 
93, 1386 (1990); M. Lozada-Cassou, W. Olivares and B. 
Sulbaran, Phys. Rev. E 53, 522 (1996). 
" M. J. Stevens and M. O. Robbins, Europhys. Lett. 12, 81 
(1990). 

I. Rouzina and V. A. Bloomfield, _J. Phys. Chem. 100 , 9977 

and J. 



(1996); see also B. I. Shklovskii, |cond-mat/9809429 



13 



J. Arenzon, J. F. Stilck and Y. Levin, |cond-mat/9806358| 
P. A. Pincus and S. A. Safran, Europhys. Lett. 42, 103 
(1998). 

^* Y. Levin, Physica A 265, 432 (1999). 
N. Ise, Angew. Chem. 25, 323 (1986). 

J. N. Israelachvih and H. K. Christenson, Physica A 140, 
278 (1986). 

J. C. Crocker and David G. Grier, Phys. Rev. Lett. 77, 
1897 (1996). 

M. D. Carbajal-Tinoco, F. Castro-Roman, and J. L. Arauz- 

Lara, Phys. Rev. E 53, 3745 (1996). 
1^ P. Tarazona, Phys. Rev. A 31, 2672 (1985). 
2" W. A. Curtin and N. W. Ashcroft, Phys. Rev. A 32, 2909 

(1985). 

A. R. Denton and N. W. Ashcroft, Phys. Rev. A 39, 4701 
(1989). 

B. Brami, J.-P. Hansen and F. Joly, Physica A 95, 505 
(1979). 

R. D. Groot, J. Chem. Phys. 95, 9191 (1991). 
2* M. J. Stevens, M. L. Falk and M. O. Robbins, J. Chem. 

Phys. 104, 5209 (1996). 
■^^ M. O. Robbins, private communications. 
2^ S. Nordholm, Chem. Phys. Lett. 105, 302 (1984). 

M. N. Tamashiro, Y. Levin and M. C. Barbosa, cond 



mat/98102l4 Physica A (in press). 
R. Penfold, S. Nordholm, B. Jonsson and C. E. Woodward, 
J. Chem. Phys. 92, 1915 (1990). 



7 



